{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notebook Examples for Chapter 6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "# these are innocuous but irritating\n",
    "warnings.filterwarnings(\"ignore\", message=\"numpy.dtype size changed\")\n",
    "warnings.filterwarnings(\"ignore\", message=\"numpy.ufunc size changed\")\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training data separability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ee\n",
    "ee.Initialize()\n",
    "\n",
    "# first 4 principal components of ASTER image\n",
    "image = ee.Image('users/mortcanty/supervisedclassification/AST_20070501_pca') \\\n",
    "            .select(0,1,2,3)\n",
    "\n",
    "# training data\n",
    "table = ee.FeatureCollection('users/mortcanty/supervisedclassification/train')\n",
    "trainData = image.sampleRegions(table,['CLASS_ID'])\n",
    "print trainData.size().getInfo()  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def jmsep(class1,class2,image,table):\n",
    "# Jeffries-Matusita separability    \n",
    "    table1 = table.filter(\n",
    "        ee.Filter.eq('CLASS_ID',str(class1-1)))\n",
    "    m1 = image.reduceRegion(ee.Reducer.mean(),table1)\\\n",
    "              .toArray() \n",
    "    s1 = image.toArray() \\\n",
    "         .reduceRegion(ee.Reducer.covariance(),table1)\\\n",
    "         .toArray()\n",
    "    table2 = table.filter(\n",
    "        ee.Filter.eq('CLASS_ID',str(class2-1)))\n",
    "    m2 = image.reduceRegion(ee.Reducer.mean(),table2)\\\n",
    "              .toArray()\n",
    "    s2 = image.toArray() \\\n",
    "        .reduceRegion(ee.Reducer.covariance(),table2,15)\\\n",
    "              .toArray()\n",
    "    m12 = m1.subtract(m2)  \n",
    "    m12 = ee.Array([m12.toList()]) # makes 2D matrix  \n",
    "    s12i = s1.add(s2).divide(2).matrixInverse()\n",
    "#  first term in Bhattacharyya distance\n",
    "    B1 = m12.matrixMultiply(\n",
    "          s12i.matrixMultiply(m12.matrixTranspose())) \\\n",
    "            .divide(8)\n",
    "    ds1 = s1.matrixDeterminant()\n",
    "    ds2 = s2.matrixDeterminant() \n",
    "    ds12 = s1.add(s2).matrixDeterminant()\n",
    "#  second term\n",
    "    B2 = ds12.divide(2).divide(ds1.multiply(ds2).sqrt())\\\n",
    "             .log().divide(2)\n",
    "    B = ee.Number(B1.add(B2).project([0]).toList().get(0))\n",
    "#  J-M separability\n",
    "    return ee.Number(1).subtract(ee.Number(1) \\\n",
    "             .divide(B.exp())) \\\n",
    "             .multiply(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print jmsep(5,9,image,table).getInfo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def band_mean(current,prev):\n",
    "    current = ee.String(current)\n",
    "    prev = ee.Dictionary(prev)\n",
    "    trainData = ee.FeatureCollection(prev.get('trainData'))\n",
    "    class_id = prev.get('class_id')\n",
    "    means = ee.List(prev.get('means'))\n",
    "    mu = trainData.filter(ee.Filter.eq('CLASS_ID',class_id)).aggregate_mean(current)\n",
    "    return ee.Dictionary({ 'trainData':trainData,'class_id':class_id,'means':means.add(mu) })\n",
    "\n",
    "def class_mean(trainData,class_id,bandNames):\n",
    "    first = ee.Dictionary({'trainData':trainData,'class_id':str(class_id),'means':ee.List([])})\n",
    "    return ee.Dictionary(bandNames.iterate(band_mean,first)).get('means')\n",
    "\n",
    "mu = ee.Array(class_mean(trainData,7,image.bandNames()))\n",
    "print mu.getInfo()    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Naive Bayes  on the GEE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import IPython.display as disp\n",
    "jet = 'black,blue,cyan,yellow,red,brown'\n",
    "\n",
    "# rename the class ids from strings to integers\n",
    "trainData = image.sampleRegions(table,['CLASS_ID'])\\\n",
    "    .remap(['0','1','2','3','4','5','6','7','8','9'],\n",
    "           [0,1,2,3,4,5,6,7,8,9],'CLASS_ID')\n",
    "    \n",
    "# train a naive Bayes classifier    \n",
    "classifier = ee.Classifier.continuousNaiveBayes()\n",
    "trained = classifier\\\n",
    "    .train(trainData,'CLASS_ID',image.bandNames())\n",
    "\n",
    "# classify the image and display    \n",
    "classified = image.classify(trained)\n",
    "url = classified.select('classification')\\\n",
    "    .getThumbURL({'min':0,'max':9,'palette':jet})\n",
    "disp.Image(url=url)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bayes Maximum Likelihood"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run scripts/classify -p [1,2,3,4] -a 1 imagery/AST_20070501_pca.tif imagery/train.shp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run scripts/dispms -f imagery/AST_20070501_pca_class.tif -c \\\n",
    "-r  \"['WATER', 'RAPESEED', 'SUGARBEET', 'SUBURBAN', 'INDUSTRIAL', 'CONIFEROUS', 'GRAIN', 'GRASSLAND', 'HERBIFEROUS', 'OPENCAST']\" \\\n",
    "-s '/home/mort/LaTeX/new projects/CRC4/Chapter6/fig6_5.eps'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gaussian kernel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run scripts/classify -p [1,2,3,4] -a 2 -P imagery/AST_20070501_pca.tif imagery/train.shp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run scripts/dispms -f imagery/AST_20070501_pca_class.tif -c \\\n",
    "-r \"['WASSER [BL', 'RAPS [YELL', 'RUEBEN [CY', 'SIEDLUNG [', 'GEWERBE [M', 'NADELWALD', 'GETREIDE [', 'GRAS [RED2', 'LAUBWALD [', 'TAGEBAU [W']\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run scripts/dispms -f imagery/AST_20070501_pca_classprobs.tif -p [4,2,1] -e 1 \\\n",
    "#-s '/home/mort/LaTeX/new projects/CRC4/Chapter6/fig6_6.eps'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ffn with backpropagation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "run scripts/classify -p [1,2,3,4] -a 3 -e 10 -L [1000] imagery/AST_20070501_pca.tif imagery/train.shp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ffn with scaled conjugate gradiant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run scripts/classify -p [1,2,3,4] -a 4 -e 1000 -L [10] imagery/AST_20070501_pca.tif imagery/train.shp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ffn with extended Kalman filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run scripts/classify -p [1,2,3,4] -a 5 -e 10 -L [10] imagery/AST_20070501_pca.tif imagery/train.shp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ffn in TensorFlow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Accessing the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import auxil.readshp as rs\n",
    "import gdal\n",
    "import matplotlib.pyplot as plt\n",
    "from osgeo.gdalconst import GA_ReadOnly\n",
    "import numpy as np\n",
    "\n",
    "# get the training data\n",
    "infile='imagery/AST_20070501_pca.tif'\n",
    "gdal.AllRegister()\n",
    "inDataset = gdal.Open(infile,GA_ReadOnly)\n",
    "pos=[1,2,3,4]\n",
    "Gs,ls,K,_ = rs.readshp('imagery/train.shp',inDataset,pos)\n",
    "m = ls.shape[0]\n",
    "print '%i observations read'%m\n",
    "\n",
    "# split into train/test\n",
    "idx = np.random.permutation(m)\n",
    "Gs = Gs[idx,:] \n",
    "ls = ls[idx,:]         \n",
    "ls = np.argmax(ls,1)\n",
    "\n",
    "Gstrn = Gs[:int(0.67*m),:]\n",
    "lstrn = ls[:int(0.67*m)] \n",
    "Gstst = Gs[int(0.67*m):,:]  \n",
    "lstst = ls[int(0.67*m):] \n",
    "\n",
    "# read entire image\n",
    "cols = inDataset.RasterXSize\n",
    "rows = inDataset.RasterYSize  \n",
    "Gs_all = np.zeros((cols*rows,4))\n",
    "k= 0\n",
    "for b in pos:\n",
    "    band = inDataset.GetRasterBand(b)\n",
    "    Gs_all[:,k] = band.ReadAsArray(0,0,cols,rows)\\\n",
    "                          .astype(float).ravel()\n",
    "    k += 1      "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Network architecture"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "from datetime import datetime\n",
    "# placeholders\n",
    "Gs = tf.placeholder(tf.float32,shape=(None,4))\n",
    "ls = tf.placeholder(tf.int64,shape=(None))\n",
    "# hidden layer with rectified linear units (relu) \n",
    "hidden=tf.layers.dense(Gs,10,activation=tf.nn.relu)\n",
    "# output layer\n",
    "logits=tf.layers.dense(hidden,10)\n",
    "# cross entropy cost function\n",
    "xentropy=tf.nn.sparse_softmax_cross_entropy_with_logits\\\n",
    "                            (labels=ls,logits=logits)\n",
    "cost=tf.reduce_mean(xentropy)\n",
    "# training algorithm with 0.01 learning rate\n",
    "optimizer=tf.train.GradientDescentOptimizer(0.01)\n",
    "training_op=optimizer.minimize(cost)\n",
    "# variables initializer \n",
    "init=tf.global_variables_initializer()\n",
    "# accuracy evaluation\n",
    "correct=tf.nn.in_top_k(logits,ls,1)\n",
    "accuracy = tf.reduce_mean(tf.cast(correct,tf.float32))\n",
    "# saver\n",
    "saver = tf.train.Saver()\n",
    "# logger for tensorboard\n",
    "cost_summary = tf.summary.scalar('COST',cost)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Training and testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "now = datetime.utcnow().strftime(\"%Y%m%d%H%M%S\")\n",
    "logdir = 'tf_logs/run-'+str(now)\n",
    "file_writer = tf.summary.FileWriter(logdir,tf.get_default_graph())\n",
    "with tf.Session() as sess:\n",
    "    init.run()\n",
    "    for epoch in range(5000):\n",
    "        if epoch % 200 ==0:\n",
    "            summary_str =cost_summary.eval(feed_dict={Gs:Gstrn,ls:lstrn})\n",
    "            file_writer.add_summary(summary_str,epoch)\n",
    "        sess.run(training_op,feed_dict={Gs:Gstrn,ls:lstrn})\n",
    "    acc = accuracy.eval(feed_dict={Gs:Gstst,ls:lstst})\n",
    "    save_path = saver.save(sess,'imagery/dnn.ckpt')\n",
    "file_writer.close()\n",
    "print 'Test accuracy: %f'%acc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!tensorboard --logdir tf_logs/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with tf.Session() as sess:\n",
    "    saver.restore(sess,'imagery/dnn.ckpt')\n",
    "    Z = logits.eval(feed_dict={Gs:Gs_all})\n",
    "    cls = np.argmax(Z,1)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 10))\n",
    "ax.imshow(np.reshape(cls/10.0,(rows,cols)),cmap='jet')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Deep learning network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run scripts/classify -p [1,2,3,4] -a 6 -e 1000 -L [10,10,10] imagery/AST_20070501_pca.tif imagery/train.shp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run scripts/ct imagery/AST_20070501_pca_Dnn(tensorflow).tst"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run scripts/dispms -f imagery/AST_20070501_pca_class.tif -c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run scripts/ct imagery/AST_20070501_pca_Dnn(tensorflow).tst"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Support Vector Machine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run scripts/classify -p [1,2,3,4] -a 7 imagery/AST_20070501_pca.tif imagery/train.shp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run scripts/ct imagery/AST_20070501_pca_SVM.tst"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# train a SVM  \n",
    "classifier = ee.Classifier.svm(kernelType='RBF',gamma=0.01,cost=100)\n",
    "trained = classifier.\\\n",
    "    train(trainData,'CLASS_ID',image.bandNames())\n",
    "\n",
    "# classify the image and display    \n",
    "classified = image.classify(trained)\n",
    "url = classified.select('classification').\\\n",
    "    getThumbURL({'min':0,'max':9,'palette':jet})\n",
    "disp.Image(url=url)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
